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^ ' Abstract 
I ' 

O ' In this paper a prescription for generating an equilibrium spherical system depicted 

H ' 

I with cuspy density profiles is extended to a core density profile, the Burkert profile, 

^ ' which is observationally more suitable to dwarfs. By using a time-saving Monte Carlo 
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method instead of N-body simulations, we show that the Burkert halo initialized with 
the distribution function that depends only on energy is in a practically stable equilib- 
rium. The one generated using the local Maxwellian approximation is unstable, where 
the fiat core density structure tends to steepen. This is fundamentally different from 
the previous study on a halo with cuspy density profile. The deviation of the unstable 
"Burkert" from the initial Burkert profile is found to be closely related to taking a 
Gaussian as the velocity distribution at any given point. The significance of not using 
the local Maxwellian approximation is further demonstrated by exploring the dynam- 
ical evolution of compact super star clusters in the Burkert halo. In particular, the 
local Maxwellian approximation will result in underestimating the dynamical friction 
and overestimating sinking time scale. Accordingly, it will lead to lower probability 
of forming massive bulges and young nuclear star clusters in bulgeless galaxies. 
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1. Introduction 

Generating initial conditions for numerical experiments of a system with a certain density 
profile is a well-defined but difficult procedure. There are two steps for constructing the initial 
conditions of one desired model: 1) calculate the steady-state distribution function of the desired 
model; 2) use Monte Carlo method to generate the initial conditions. The main difficulty comes 
from the first step. One possibility is to make some simplifying approximations about the nature 
of the distribution function (Toomre, & Toomre 1972; Hernquist & Quinn 1988). Hernquist 
(1993) described a prescription for constructing N-body realization assuming that velocity 
distribution at any given point is Maxwellian. It is called the local Maxwellian approximation 
(Kazantzidis, Magorrian & Moore 2004), of which the advantage is easy to implement. However, 
Kazantzidis, Magorrian & Moore (2004, KMM hereafter) show that the halo initialized using 
the local Maxwellian approximation (the Maxwellian halo hereafter for simplicity ) can be 
significantly far from equilibrium, and that for high resolution simulations it is more reasonable 
to use a stable system generated with a certain steady-state distribution function (the DF halo 
hereafter) . 

In particular, they demonstrated that an unstable cuspy halo constructed using the local 
Maxwellian approximation will soon lose its cusp structure due to its own evolution, which is 
previously claimed to be due to minor merger between a host halo and the falling sub-system. 
On the contrary, a stable cuspy halo survives the minor merger. This investigation indicates 
that the non-equilibrium Maxwellian halo will surely lead to spurious evolution of the host 
halo when dealing with the response, e.g. the responding of a host halo to sinking SSCs will 
otherwise be mixed with the artificial evolution of the host halo, induced by its unstability. 
This effect has been investigated in our recent work on the bulge formation in late-type spirals 
(Fu et al 2005). 

The halos on which KMM studied are those with cuspy structure. They have not 
explored the halos with flat inner profile, e.g. the Burkert profile. The progress in observations 
(e.g. Marchesini et al. 2002; and the references therein) shows, however, that the cuspy DM 
profiles, such as the NFW (Navarro, Frenk & White 1997) profile, are not adequate for a large 
fraction of dwarf galaxies, which are dominated by DM halos, and suggests that the core density 
profiles, the Burkert (1995) profile for example, are more suitable to these galaxies. 

It would then be very much instructive to extend KMM's study to the Burkert halo. It 
is the halo that contains less mass in the central region than those in the cuspy halos, e.g. the 
NFW profile. It follows that the technical difficulties will arise for high-resolution simulations 
of the core halos, especially for those studies focused on the central regions (Fu, Huang & 



2 



Deng 2003a; Huang, Deng & Fu 2003). For a given mass of DM halo, our test shows that in 
order to reach the same resolution in the central lOpc region the number of particle for the 
Burkert halo should be a hundred times more than that for the NFW halo. It would need much 
more computing time with the N-body simulations. 

In this paper, we aim at describing the way to generate the initial conditions for an 
isotropic, spherical Burkert halo, and perform simulations using the Monte Carlo method 
(Henon 1971a,b) to study the stable situation of the Maxwellian and the DF Burkert ha- 
los. This may have important implications for the dynamical evolution of SSCs in a responding 
dark matter halo, resulting in different bulge formation history in late-type spirals. 

The organization of the paper is as follows. The procedure to constructing the initial 
conditions for our model is described in section 2. Section 3 gives the results of the simulations. 
We discuss the results and illustrate the important implications of generating the DF halo in 
Section 4. Finally, Section 5 notes our conclusion. 

2. Models and Methods 

2.1. Density Profile 

Here we study the Burkert dark matter halo (1995): 

P (^i + r/rs)[l + ir/r,y] 
where the central density ps is taken to be O.OSMq/^c^ following Marchesini et al. (2002), and 
the scale radius, is computed from the mass of the DM halo M200 = IO^^Mq. In the inner part 
of the Bukert halo the profile has a core structure, while the slope approximates to —3 in the 
infinity. The cumulative mass profile with such a distribution diverges as r — 00. In fact, the 
Burkert profile is a phenomenological formula, which provides a good fit to the observed data 
to less than the virial radius r^r and not valid to arbitrarily large distance from the galactic 
center. Instead of truncating the profile sharply at r^r, we have chosen an exponential cut-off 
for r > following KMM, which sets in at r^r and turn-off in a scale r decay 

P= ,.^ ,A —rexp[-'-^] for{r>r^,r) (2) 

\i- I Cjyl -r C J Tuij, f decay 

where c = r^ir/rs is the concentration parameter. In order to ensure a smooth transition between 
eq.(l) and eq.(2) at r^ir, we require the logarithmic slope there to be continuous. This leads to 

decay 1 ~t~ C 1 ~1~ C 

We consider the simulations between a minimum radius Tmin and a maximum radius Vmax- The 
minimum radius is chosen such that it has sufficient particles for us to analyze; The maximum 
radius is equal to the virial radius plus several r decay Here we adopt three r decay i^decay = O.lr^j^). 
In addition, we add a boundary at the maximum radius. If the halo is stable, it is in homeostasis, 
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that is, at any time the number of particles which move out off the boundary is equal to that 
of particles which move in. In the case of the Burkert halo, the DF (Burkert) halo will not be 
stable in the outer region, if not adding such a boundary. This is because, not like the case of 
KMM, the mass outside Tmr+srdecay is still not negligible. In other words, a sharp truncation 
there will make the DF halo unstable. 

2.2. Distribution Function 

In this paper, we restrict our models to isotropic, non-rotating halos. According to the 
Jeans' theorem (Lynden-Bell 1962; Binney & Tremaine 1987), the distribution function of any 
steady-state spherical system can be expressed as f{E,L), where E is the binding energy and 
L is the angular momentum vector. As usual, the relative potential ^' and the relative energy 
e (a dimensionless energy in units of n'^Gpsr'^) are defined as: 



where $ is the potential and $0 is a constant fulfilling / > for e > and / = for e < 0. 
In an isotropic spherical model the distribution function depends on e only, that is, f = f{e). 
Integrating / over all velocities, we can get the density profile: 



The inversion of the above equation gives the distribution function (Eddington 1916; Binney 
& Tremaine 1987), 



The (Pp/d^"^ factor can be evaluated from eq.(l) and eq.(2). The second term of the right-hand 
side vanishes in large distance in our model. 

2.3. Initialization Procedure 

Here we describe two approaches to initialize the numerical simulations for the Burkert 
profile using the local Maxwellian approximation and the steady-state distribution function, 
respectively. 

2.3.1. Exact Distribution Function 

Having the density profile p(r), we can calculate the model's cumulative mass distribu- 
tion M(r) and the gravitational potential $(r). Then we can calculate the distribution function 
from eq.(6). 

The integrand in eq.(6) diverges at one or both of the limits, but this can be solved 
using standard techniques for improper integral (Press et al 1986). The distribution function 
is obtained on a grid of e equally spaced in log(£). The accuracy of the numerical integration 
is checked by comparing the density profile derived from eq.(5) and that given by eq.(l). The 
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two profiles are very well accordant. We find that the distribution function is non- negative 
everywhere, which proves that the assumption of the isotropic velocity distribution for the 
spherical Burkert model is reasonable and physical. 

The procedure to obtain the physical quantities of the sampling particles is the following. 
We randomly sample the initial positions and velocities of the particles. For our model, the 
probability density function (Windrow 2000) of the particle having the relative energy e at the 
radius R is: 

P{e,R)^R\^>{R)-efl'f{e) (7) 

We use the rejection method (Press et al 1986) to sample the particles in two steps: 1) 
Generate the position of the particle. The probability density function with radius R is: 

P{R) oc R^p{R) (8) 

2) Generate the energy of the particle. Once given R, we can get the relative potential '^{R)- 
Then the probability density function having the relative energy e and potential '^{R) is: 

P{e)^{^{R)-ef''f{e) (9) 

Once the position R and the relative energy e of a particle are given, the speed v can be easily 
determined by e=v'^ 12+^^ . Finally, the tangential and radial velocities of each particle can be 
randomly sampled considering that our model is isotropic. 

2.3.2. Local Maxwellian Approximation 

The way to sample the position of a particle is the same as that depicted above. The 
difference is how to generate the velocity of the particle. The collisionless Boltzmann equation 
for a spherical system can be written as: 

^ + f|2^?-W + ^|^-,f (10) 

where v"^, Vg, and v"^ are the velocity dispersions in spherical coordinates (Binney & Tremaine 
1987). Assuming the model is isotropic, it can be integrated to: 
1 roo 

^r = ^^l P(r)^dr (11) 

Using the known density profile, the radial dispersion can be computed as a function of ra- 
dius. In the local Maxwellian approximation, the ID velocity distribution at a given point is 
approximated by a Gaussian (Hernquist 1993): 

F{v,r) = 47T{-^)^/\'^exp{-vy2^) (12) 

The speed of a halo particle is initialized from the eq.(12). Then the tangential and radial 
velocities are randomly sampled in the same way as that described in the above subsection. 
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2.4- Simulations and Parameters 

Many approaches have been used in simulations to study the galactic dynamics. One is to 
give the statistical description of the system, e.g. distribution function f{1^,lf,t). Rosenbluth 
et al. (1957) discussed a way to calculate / by solving a much complicated equation, Fokker- 
Planck equation. But it needs to make a number of arbitrary simplifications. A commonly used 
alternative way is N-body simulations. Unfortunately, this approach requires much computing 
time, which will greatly increase with the number of particles. As we mentioned before, in the 
case of a core halo like the Burkert one (1995), a large number of particles are in demand to 
keep high resolution in the central region. In view of this, we adopt a time-saving method, the 
Monte Carlo procedure (Henon 1971a,b), rather than the N-body simulations. This procedure 
has been proved to reproduce the behavior of the system given by the Fokker-Planck equation 
and to be much faster than the N-body simulations. 

The basic ideas of the Monte Carlo procedure is the following. We can divide the 
gravitation filed of the system into two parts: a main smoothed-out field and a small fluctuating 
one. The particle moves during every time step At which satisfies: 

T, < At < Tr (13) 

where Tc and Tr are the crossing time and the relaxation time, respectively. The halo evolves 
accordingly. The relaxation time is related to the crossing time (Binney & Tremaine 1987): 

Tr/Tc(x N/ In N (14) 

And the crossing time can be calculated as a function of the total mass M and the total energy 
e of the system (van Albada 1968): 

T, = CGM''/^\e\-^/^ (15) 

where C is a numerical constant. 

Neglecting the fluctuating field in a first approximation, the motion of the particle is 
governed by the main spherical symmetry field which changes with the density or mass distribu- 
tion of the halo. The particle will then take a plane rosette motion (Binney & Tremaine 1987). 
In any given state the particle is characterized by its position, energy and its angular momen- 
tum. One could derive the statistic of the system by randomly choosing the new positions and 
velocities of the particles on their respective orbits. Taking them as the new conditions, we can 
calculate a new density (or mass) distribution and accordingly a new potential of the system, 
which will in turn determine the positions of particles at the next step. The equilibrium state 
of the halo can be tested by comparing the new density (or mass) distribution with the initial 
one. Due to this character, if we find some deviation from the initial state of a system at one 
step, the deviation should grow in the simulations until the system relaxes to a new stable 
equilibrium (as Figure lb shows). In other words, what we could illustrate with the Monte 
Carlo method at the present time is that the system is practically stable or unstable. 
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Table 1. Parameters of Simulations 
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Fig. 1. Radial density profiles for Model A and B are presented in a) and b) panel, respectively. The 
dotted, dashed and dot-dashed lines in both panels correspond to the density distributions derived at 
three randomly choosing steps. The thin solid lines in two panels show the initial Burkert profile given by 
eq.(l). The insert diagram in b) illustrates the situation in the central region for clarity. 

For our simulations, the relaxation time of the system is much long and we need not 
consider the two-body relaxation in our study. If the system is in stable equilibrium, the 
statistic of the system such as the potential field (or the density and mass distribution) will 
keep the same though the positions and velocities of the particles are randomly chosen in each 
step. On the contrary, in an unstable system, from the density or mass profile we can know 
that the statistic of the system have changed after one step. 

The parameters of our simulations are given in Table 1. 

3. Results 



We run our models for 20 steps and generate 20 individual snapshots. The radial density 
profiles shown in this section are obtained by binning the particles from individual snapshot. 
Figure la) illustrate the radial density profiles for Model A, where the density distributions 
at three different steps, denoted with dotted, dashed and dot-dashed lines, are plotted. For 



7 



comparison, the (initial) Burkert profile given by eq.(l) is also shown in the figure, denoted 
with thin solid line. Obviously the radial density profiles for Model A have no evolution over 
the simulations, all of them are indeed very well matched to the (initial) Burkert profile. In 
view of this, we conclude that the DF Burkert halo, i.e. the Burkert halo initialized with the 
exact distribution function, is in stable equilibrium. 

The opposite situation occurs for Model B. We present the results in Figure lb), where 
the density distributions derived at three randomly choosing steps are indicated. The profiles 
denoted with dotted, dashed and dot-dashed lines illustrate significant deviation from the (ini- 
tial) Burkert profile at the central and outer regions. The central distribution tends to steepen 
as clearly shown in the insert diagram of Figure lb) and the outer densities decline. The dis- 
parity from the initial one simply shows that the Maxwellian Burkert halo, i.e. the Burkert 
halo generated with the local Maxwellian approximation is not in equilibrium. 

4. Discussions 

Following KMM, we analyze the actual velocity structure at various distances for Model 
A, compared with the corresponding Gaussian distributions used in the Maxwellian approxi- 
mation (Model B). The solid and dotted lines in Figure 2 correspond to the true and Gaussian 
velocity distributions, respectively, at four different distances from the center. It is evident that 
only at about the scale radius the Gaussian distribution is a close approximation to the reality. 
The disparity in the velocity structure from the Gaussian distribution is very strong near the 
center. The true velocity distributions there are strongly peaked than a Gaussian. It is still 
apparent at far from the center, where the true velocity structure is shallower than a Gaussian 
on the contrary. 

The feature in velocity distributions described above is indeed refiected in the Maxwellian 
Burkert halo, where the deviation from the (initial) Burkert halo is obvious at the central and 
outer regions. The calculations show that the total energy and angular momentum of particles 
in the Maxwellian Burkert halo are less than those in the DF Burkert halo, i.e. the (initial) 
Burkert halo. That is, the Maxwellian Burkert halo is colder than the DF Burkert halo. This 
explains the hoist of the central density profile in the Maxwellian halo. In other words, the 
unstable situation of the Maxwellian halo is closely related to taking a Gaussian as the velocity 
distribution at any given point. 

In some existing N-body simulations, the considered system with a predetermined mass 
density profile are initialized using the Maxwellian approximation. In order to start simulations 
with a system in stable or quasi-stable equilibrium, the authors allow the constructed system to 
relax for some time so that a ceratin quasi-stable equilibrium is attained. According to KMM's 
investigation and the study presented in this work, however, the simulations should start with a 
system not having the predetermined mass density profile. That is to say, spurious "evolution" 
of the system initialized with the Maxwellian approximation still cannot be avoided in their 
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Fig. 2. Histograms of the one dimensional velocity distribution at four differ- 
ent distances from the center in units of kpc. The solid lines correspond to the 
true velocity distribution, and the dotted lines show the Gaussian velocity profile. 

way. 

Here we would like to further demonstrate the significance of initializing a system with 
the distribution function by exploring the dynamical evolution of a compact super star cluster 
(SSC hereafter) in a DM dominated galaxy, depicted with the Burkert profile. The term 
"compact" here means no tidal stripping is effected. In this case, the most important process 
involved is the dynamical friction the compact SSC experiences, which reads as (Binney & 
Tremaine 1987) 



dV 



M 



dt 



-16n^\nAG^m{M + m 



■V 



M 



where 



A 



h 

'-'max ' typ 



(16) 



(17) 



G{M + m) 

M and Vm (with Vm = \^m\) are, respectively, the mass and velocity of the SSC experiencing 
the dynamical friction. The mass of the DM particle, m, is much smaller than M. The quantity 
bmax is the maximum impact parameter and Vtyp a typical speed. Neither bmax nor Vtyp is 
precisely defined. Following Binney & Tremaine (1987), we take bmax = 2 kpc and Vtyp = Vm- 

If the DM has a Maxwellian velocity distribution with dispersion abkgd, the dynamical 
friction can be written as: 
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dVM 27rlog(l + A2)G2 



Mp 



[erf{X) 



2X 



exp(-X2)]VM 



(18) 



where erf is the error function, and 



(19) 




The velocity dispersion ahkgd can be evaluated from the Jeans equation. 

Eq.(16) and eq.(18) are, respectively, the basis for considering the dynamical friction in 
the DF Burkert halo and the Maxwellian Burkert one. The initial SSC is assumed to move at 
a radius r=lkpc with the local circular velocity. The sinking history of the SSC is presented in 
Figure 3 a) and b) for SSC with different masses. The thick solid lines in both panels express 
the case in Model A, i.e. in the Burkert halo initialized using the distribution function, while 
the thick dotted lines the case in Model B, i.e. in the Maxwellian Burkert halo. 

We can see from Figure 3a) that the SSC with typical mass of 2 x IO^Mq (Fu, Huang 
& Deng 2003a) sinks quickly, in about 4 x lO^yr, to the central region within a hundred pc 
in Model A, while takes much longer time in Model B. It means that the SSC in Model A 
experiences much stronger dynamical friction and loses angular momentum quickly at first. 
This situation is consistent with the velocity distributions that we discussed above. In fact, 
in an isotropic halo, the effect of the dynamical friction comes from the DM particles whose 
velocities are less than that of the SSC (Binney & Tremaine 1987). As we illustrated in Figure 2, 
the true velocity structure in the inner part of the halo is more peaked than a Gaussian. Thus 
more particles with low velocities exist, resulting in stronger dynamical friction on SSC. Inside 
200 pc from the center, the SSC in Model A sinks slowly than that in Model B, mainly due to 
its low velocity and low deceleration accordingly. 

A similar sinking story is presented in Figure 3b) for a SSC with mass of 8 x IO^Mq, 
where we can find same trend of faster sinking in Model A than in Model B. The interesting 
thing is that the SSC with smaller mass reaches the center, less than 10 pc, in shorter period 
of time, as compared with the corresponding situation in Model A shown in Figure 3a). The 
compact SSC with larger mass experiences stronger dynamical friction, resulting in stronger 
deceleration but longer sinking time. 

As we know, when the dynamical friction exists the periodic motion of a test particle 
in a spherical potential is damped like that of a damped simple harmonic oscillator. If the 
dynamical friction is strong enough, it is possible for the motion to lose completely its oscillation 
behavior, the so-called over-damping. For an over-damped SSC, the energy of SSC is in the 
form of potential energy in most of time, which makes the energy loss by dynamical friction 
much less effective than otherwise. For an incompact SSC, the tidal stripping must not be 
ignored. In that case, the peeling process will make a heavy SSC lighter, and the over-damping 
may not appear. 

The shorter sinking time for SSCs in a DM halo initialized with the steady-state distri- 
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Fig. 3. Sinking history of a compact SSC in the Burkert halo. Panel a) and b) present the situation for 
SSC with mass of 2 x 10^ Mq and 8 x 10^ Mq, respectively. The thick solid lines in both panels correspond 
to the case in Model A, i.e. the Burkert halo initialized with the distribution function, while the thick 
dotted ones in Model B: the Maxwellian halo. The thin solid line denotes the distance of 10 pc to the 
center 

bution function should lead to higher probability of forming massive bulges in late-type galaxies 
(Fu, Huang & Deng 2003a), as well as of forming nuclear star clusters at young ages (Huang, 
Deng & Fu 2003). It will also bring the intermediate mass black holes, the seed black holes, 
within SSCs to the center of a DM halo at early stage (Fu, Huang & Deng 2003b). All of these 
matter are hot topics in astrophysics, which emphasizes the important implications of using 
the distribution function to initialize the system. 

5. CONCLUSION 

Constructing appropriate and reasonable initial conditions for an isolated equilibrium 
system is one of the most important points to various numerical simulations related to the 
formation and evolution of galaxies. In this paper, we extend KMM's study on cuspy density 
profiles to the Burkert profile. Using a time-saving Monte Carlo method, we have shown clearly 
that the Burkert halo initialized with the local Maxwellian approximation tends to steepen, 
which is closely related to adopting a Gaussian as the velocity distribution at any given points. 
This spurious evolution leads to underestimating the dynamical friction on SSCs moving in dark 
matter dominated galaxies. This important demonstration gives us a valuable clue to clarify 
the embarrassed situation in our previous investigation (Fu, Huang & Deng 2003a; Huang, 
Deng & Fu 2003), i.e. the adopted Maxwellian approximation in these works causes longer 
sinking time for SSCs so as not to form massive bulges in a few Gyrs and not to form a young 
nuclear star cluster in the Burkert halo. To initialize the dark matter halos with the steady- 
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state distribution function should be critical to what we are pursuing, as what we have shown 
in this paper. 

The authors would like to thank the referee, Audi Burkert, for his valuable comments which 
improve the description of the paper. This work is supported by NKBRSF G19990754 and 
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